#Roland Neil and Zach Wehrwein
#Gov 2001 Paper
#This file carries out the analysis for the final paper.      
#April 19th, 2016

setwd("C:/Users/Roland/Documents/Harvard/Gov2001/Paper")
rm(list=ls()) #shortcut to remove data

library(readstata13)
library(lmtest)
library(sandwich)
library(dplyr)

######Table 1#####
all_atus <- read.dta13("New Data/all_atus.dta")
all_atus<-all_atus[!(all_atus$age<18 | all_atus$age>65 | all_atus$unclassified_paper>0),]

#Column 1 (same as table 1, column 3 in AER paper)
out1 <- lm(work_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out1, vcov = vcovHC(out1, "HC1")) 
out2 <- lm(worka_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out2, vcov = vcovHC(out2, "HC1")) 
out3 <- lm(worku_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out3, vcov = vcovHC(out3, "HC1")) 
out4 <- lm(childcare_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out4, vcov = vcovHC(out4, "HC1")) 
out5 <- lm(home_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out5, vcov = vcovHC(out5, "HC1")) 
out6 <- lm(homeproduction_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out6, vcov = vcovHC(out6, "HC1")) 
out7 <- lm(homeown_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out7, vcov = vcovHC(out7, "HC1")) 
out8 <- lm(shopping_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out8, vcov = vcovHC(out8, "HC1")) 
out9 <- lm(othercare_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out9, vcov = vcovHC(out9, "HC1")) 
out10 <- lm(leisure_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out10, vcov = vcovHC(out10, "HC1")) 
out11 <- lm(tv_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out11, vcov = vcovHC(out11, "HC1")) 
out12 <- lm(socializing_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out12, vcov = vcovHC(out12, "HC1")) 
out13 <- lm(sleeping_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out13, vcov = vcovHC(out13, "HC1")) 
out14 <- lm(ep_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out14, vcov = vcovHC(out14, "HC1")) 
out15 <- lm(otherleisure_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out15, vcov = vcovHC(out15, "HC1")) 
out16 <- lm(other_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out16, vcov = vcovHC(out16, "HC1")) 
out17 <- lm(education_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out17, vcov = vcovHC(out17, "HC1")) 
out18 <- lm(civic_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out18, vcov = vcovHC(out18, "HC1")) 
out19 <- lm(ownmedical_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
coeftest(out19, vcov = vcovHC(out19, "HC1")) 

#Column 2 (average time use for 2011/2012)
out1 <- lm(work_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out1, vcov = vcovHC(out1, "HC1")) 
out2 <- lm(worka_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out2, vcov = vcovHC(out2, "HC1")) 
out3 <- lm(worku_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out3, vcov = vcovHC(out3, "HC1")) 
out4 <- lm(childcare_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out4, vcov = vcovHC(out4, "HC1")) 
out5 <- lm(home_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out5, vcov = vcovHC(out5, "HC1")) 
out6 <- lm(homeproduction_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out6, vcov = vcovHC(out6, "HC1")) 
out7 <- lm(homeown_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out7, vcov = vcovHC(out7, "HC1")) 
out8 <- lm(shopping_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out8, vcov = vcovHC(out8, "HC1")) 
out9 <- lm(othercare_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out9, vcov = vcovHC(out9, "HC1")) 
out10 <- lm(leisure_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out10, vcov = vcovHC(out10, "HC1")) 
out11 <- lm(tv_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out11, vcov = vcovHC(out11, "HC1")) 
out12 <- lm(socializing_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out12, vcov = vcovHC(out12, "HC1")) 
out13 <- lm(sleeping_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out13, vcov = vcovHC(out13, "HC1")) 
out14 <- lm(ep_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out14, vcov = vcovHC(out14, "HC1")) 
out15 <- lm(otherleisure_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out15, vcov = vcovHC(out15, "HC1")) 
out16 <- lm(other_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out16, vcov = vcovHC(out16, "HC1")) 
out17 <- lm(education_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out17, vcov = vcovHC(out17, "HC1")) 
out18 <- lm(civic_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out18, vcov = vcovHC(out18, "HC1")) 
out19 <- lm(ownmedical_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
coeftest(out19, vcov = vcovHC(out19, "HC1")) 


#Column 3 (average time use for 2013/2014)
out1 <- lm(work_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out1, vcov = vcovHC(out1, "HC1")) 
out2 <- lm(worka_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out2, vcov = vcovHC(out2, "HC1")) 
out3 <- lm(worku_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out3, vcov = vcovHC(out3, "HC1")) 
out4 <- lm(childcare_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out4, vcov = vcovHC(out4, "HC1")) 
out5 <- lm(home_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out5, vcov = vcovHC(out5, "HC1")) 
out6 <- lm(homeproduction_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out6, vcov = vcovHC(out6, "HC1")) 
out7 <- lm(homeown_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out7, vcov = vcovHC(out7, "HC1")) 
out8 <- lm(shopping_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out8, vcov = vcovHC(out8, "HC1")) 
out9 <- lm(othercare_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out9, vcov = vcovHC(out9, "HC1")) 
out10 <- lm(leisure_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out10, vcov = vcovHC(out10, "HC1")) 
out11 <- lm(tv_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out11, vcov = vcovHC(out11, "HC1")) 
out12 <- lm(socializing_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out12, vcov = vcovHC(out12, "HC1")) 
out13 <- lm(sleeping_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out13, vcov = vcovHC(out13, "HC1")) 
out14 <- lm(ep_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out14, vcov = vcovHC(out14, "HC1")) 
out15 <- lm(otherleisure_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out15, vcov = vcovHC(out15, "HC1")) 
out16 <- lm(other_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out16, vcov = vcovHC(out16, "HC1")) 
out17 <- lm(education_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out17, vcov = vcovHC(out17, "HC1")) 
out18 <- lm(civic_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out18, vcov = vcovHC(out18, "HC1")) 
out19 <- lm(ownmedical_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
coeftest(out19, vcov = vcovHC(out19, "HC1")) 


#Column 4: Comparing (average) weekly time use in 09/10 to 13/14
all_atus$dummyrecovery <- ifelse(all_atus$year<=2012, 0, 1)
out1 <- lm(work_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out1, vcov = vcovHC(out1, "HC1")) 
out2 <- lm(worka_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out2, vcov = vcovHC(out2, "HC1")) 
out3 <- lm(worku_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out3, vcov = vcovHC(out3, "HC1")) 
out4 <- lm(childcare_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out4, vcov = vcovHC(out4, "HC1")) 
out5 <- lm(home_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out5, vcov = vcovHC(out5, "HC1")) 
out6 <- lm(homeproduction_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out6, vcov = vcovHC(out6, "HC1")) 
out7 <- lm(homeown_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out7, vcov = vcovHC(out7, "HC1")) 
out8 <- lm(shopping_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out8, vcov = vcovHC(out8, "HC1")) 
out9 <- lm(othercare_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out9, vcov = vcovHC(out9, "HC1")) 
out10 <- lm(leisure_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out10, vcov = vcovHC(out10, "HC1")) 
out11 <- lm(tv_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out11, vcov = vcovHC(out11, "HC1")) 
out12 <- lm(socializing_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out12, vcov = vcovHC(out12, "HC1")) 
out13 <- lm(sleeping_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out13, vcov = vcovHC(out13, "HC1")) 
out14 <- lm(ep_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out14, vcov = vcovHC(out14, "HC1")) 
out15 <- lm(otherleisure_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out15, vcov = vcovHC(out15, "HC1")) 
out16 <- lm(other_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out16, vcov = vcovHC(out16, "HC1")) 
out17 <- lm(education_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out17, vcov = vcovHC(out17, "HC1")) 
out18 <- lm(civic_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out18, vcov = vcovHC(out18, "HC1")) 
out19 <- lm(ownmedical_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out19, vcov = vcovHC(out19, "HC1")) 

#Column 5: same as 4, with controls (conditional model)
all_atus$age1 <- 0
all_atus$age1[all_atus$age >=18 & all_atus$age <=27] <- 1
all_atus$age2 <- 0
all_atus$age2[all_atus$age >=28 & all_atus$age <=37] <- 1
all_atus$age3 <- 0
all_atus$age3[all_atus$age >=38 & all_atus$age <=47] <- 1
all_atus$age4 <- 0
all_atus$age4[all_atus$age >=48 & all_atus$age <=57] <- 1
all_atus$age5 <- 0
all_atus$age5[all_atus$age >=58 & all_atus$age <=65] <- 1
all_atus$educ1 <- 0
all_atus$educ1[all_atus$grade<12] <- 1
all_atus$educ2 <- 0
all_atus$educ2[all_atus$grade==12] <- 1
all_atus$educ3 <- 0
all_atus$educ3[all_atus$grade>12 & all_atus$grade<16] <- 1
all_atus$educ4 <- 0
all_atus$educ4[all_atus$grade>=16] <- 1

out1 <- lm(work_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out1, vcov = vcovHC(out1, "HC1")) 
out2 <- lm(worka_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out2, vcov = vcovHC(out2, "HC1")) 
out3 <- lm(worku_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out3, vcov = vcovHC(out3, "HC1")) 
out4 <- lm(childcare_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out4, vcov = vcovHC(out4, "HC1")) 
out5 <- lm(home_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out5, vcov = vcovHC(out5, "HC1")) 
out6 <- lm(homeproduction_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out6, vcov = vcovHC(out6, "HC1")) 
out7 <- lm(homeown_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out7, vcov = vcovHC(out7, "HC1")) 
out8 <- lm(shopping_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out8, vcov = vcovHC(out8, "HC1")) 
out9 <- lm(othercare_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out9, vcov = vcovHC(out9, "HC1")) 
out10 <- lm(leisure_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out10, vcov = vcovHC(out10, "HC1")) 
out11 <- lm(tv_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out11, vcov = vcovHC(out11, "HC1")) 
out12 <- lm(socializing_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out12, vcov = vcovHC(out12, "HC1")) 
out13 <- lm(sleeping_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out13, vcov = vcovHC(out13, "HC1")) 
out14 <- lm(ep_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out14, vcov = vcovHC(out14, "HC1")) 
out15 <- lm(otherleisure_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out15, vcov = vcovHC(out15, "HC1")) 
out16 <- lm(other_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out16, vcov = vcovHC(out16, "HC1")) 
out17 <- lm(education_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out17, vcov = vcovHC(out17, "HC1")) 
out18 <- lm(civic_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out18, vcov = vcovHC(out18, "HC1")) 
out19 <- lm(ownmedical_paper ~ dummyrecovery + male + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
coeftest(out19, vcov = vcovHC(out19, "HC1")) 

#Column 6: difference between trend (03-08) and 13/14 observation (i.e. column 3)
#Need trend values to complete, then substract those from column 3 values

######Table 2#######
#Do not run stuff above after this is run, as I am changing the main dataset (excluding males)
#This code is identical to that above, expect column 5 which does not control for "male" since all are male
#I have excluded robust standard error code for brevity, it is the same as above if needed
all_atus <- all_atus[all_atus$male==1,]

#Column 1 (average time use for 2009/2010)
out1 <- lm(work_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out2 <- lm(worka_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out3 <- lm(worku_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out4 <- lm(childcare_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out5 <- lm(home_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out6 <- lm(homeproduction_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out7 <- lm(homeown_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out8 <- lm(shopping_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out9 <- lm(othercare_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out10 <- lm(leisure_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out11 <- lm(tv_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out12 <- lm(socializing_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out13 <- lm(sleeping_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out14 <- lm(ep_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out15 <- lm(otherleisure_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out16 <- lm(other_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out17 <- lm(education_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out18 <- lm(civic_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)
out19 <- lm(ownmedical_paper ~ 1, data=subset(all_atus, year>=2009 & year<=2010), weights=weight_adj)

#Column 2 (average time use for 2011/2012)
out1 <- lm(work_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out2 <- lm(worka_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out3 <- lm(worku_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out4 <- lm(childcare_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out5 <- lm(home_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out6 <- lm(homeproduction_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out7 <- lm(homeown_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out8 <- lm(shopping_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out9 <- lm(othercare_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out10 <- lm(leisure_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out11 <- lm(tv_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out12 <- lm(socializing_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out13 <- lm(sleeping_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out14 <- lm(ep_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out15 <- lm(otherleisure_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out16 <- lm(other_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out17 <- lm(education_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out18 <- lm(civic_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)
out19 <- lm(ownmedical_paper ~ 1, data=subset(all_atus, year>=2011 & year<=2012), weights=weight_adj)


#Column 3 (average time use for 2013/2014)
out1 <- lm(work_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out2 <- lm(worka_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out3 <- lm(worku_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out4 <- lm(childcare_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out5 <- lm(home_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out6 <- lm(homeproduction_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out7 <- lm(homeown_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out8 <- lm(shopping_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out9 <- lm(othercare_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out10 <- lm(leisure_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out11 <- lm(tv_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out12 <- lm(socializing_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out13 <- lm(sleeping_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out14 <- lm(ep_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out15 <- lm(otherleisure_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out16 <- lm(other_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out17 <- lm(education_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out18 <- lm(civic_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)
out19 <- lm(ownmedical_paper ~ 1, data=subset(all_atus, year>=2013 & year<=2014), weights=weight_adj)

#Column 4
out1 <- lm(work_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out2 <- lm(worka_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out3 <- lm(worku_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out4 <- lm(childcare_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out5 <- lm(home_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out6 <- lm(homeproduction_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out7 <- lm(homeown_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out8 <- lm(shopping_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out9 <- lm(othercare_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out10 <- lm(leisure_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out11 <- lm(tv_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out12 <- lm(socializing_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out13 <- lm(sleeping_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out14 <- lm(ep_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out15 <- lm(otherleisure_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out16 <- lm(other_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out17 <- lm(education_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out18 <- lm(civic_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)
out19 <- lm(ownmedical_paper ~ dummyrecession, data=subset(all_atus, year>=2006), weights=weight_adj)

#Column 4: Comparing (average) weekly time use in 09/10 to 13/14
out1 <- lm(work_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out2 <- lm(worka_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out3 <- lm(worku_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out4 <- lm(childcare_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out5 <- lm(home_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out6 <- lm(homeproduction_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out7 <- lm(homeown_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out8 <- lm(shopping_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out9 <- lm(othercare_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out10 <- lm(leisure_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out11 <- lm(tv_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out12 <- lm(socializing_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out13 <- lm(sleeping_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out14 <- lm(ep_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out15 <- lm(otherleisure_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out16 <- lm(other_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out17 <- lm(education_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out18 <- lm(civic_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out19 <- lm(ownmedical_paper ~ dummyrecovery, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)

#Column 5: same as 4, with controls (conditional model)
out1 <- lm(work_paper ~ dummyrecovery + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out2 <- lm(worka_paper ~ dummyrecovery + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out3 <- lm(worku_paper ~ dummyrecovery  + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out4 <- lm(childcare_paper ~ dummyrecovery + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out5 <- lm(home_paper ~ dummyrecovery  + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out6 <- lm(homeproduction_paper ~ dummyrecovery + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out7 <- lm(homeown_paper ~ dummyrecovery  + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out8 <- lm(shopping_paper ~ dummyrecovery  + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out9 <- lm(othercare_paper ~ dummyrecovery + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out10 <- lm(leisure_paper ~ dummyrecovery  + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out11 <- lm(tv_paper ~ dummyrecovery  + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out12 <- lm(socializing_paper ~ dummyrecovery + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out13 <- lm(sleeping_paper ~ dummyrecovery  + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out14 <- lm(ep_paper ~ dummyrecovery + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out15 <- lm(otherleisure_paper ~ dummyrecovery + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out16 <- lm(other_paper ~ dummyrecovery + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out17 <- lm(education_paper ~ dummyrecovery + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out18 <- lm(civic_paper ~ dummyrecovery + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)
out19 <- lm(ownmedical_paper ~ dummyrecovery + black + married + hv_child + age1 + age2 + age3 + age4 + age5 + educ1 + educ2 + educ3 + educ4, data=subset(all_atus, year==2009| year== 2010| year==2013|year==2014), weights=weight_adj)

#Column 6: difference between trend (03-08) and 13/14 observation (i.e. column 3)
#Need trend values to complete, then substract those from column 3 values.

######Table 3######
#Note that all numbers in the AER article are multiplied by 100.
#Table 3, column 1: 
#First, I need appendix B1 column 1:
all_atus <- read.dta13("New Data/all_atus.dta")
all_atus<-all_atus[(all_atus$age>=18 & all_atus$age<=65 & all_atus$other==0),]
a<-weighted.mean(all_atus$work_paper, w=all_atus$weight_adj)
b<-weighted.mean(all_atus$worka_paper, w=all_atus$weight_adj)
c<-weighted.mean(all_atus$worku_paper, w=all_atus$weight_adj) 
d<-weighted.mean(all_atus$childcare_paper, w=all_atus$weight_adj)
e<-weighted.mean(all_atus$home_paper, w=all_atus$weight_adj) 
f<-weighted.mean(all_atus$homeproduction_paper, w=all_atus$weight_adj)
g<-weighted.mean(all_atus$homeown_paper, w=all_atus$weight_adj)  
h<-weighted.mean(all_atus$shopping_paper, w=all_atus$weight_adj)   
i<-weighted.mean(all_atus$othercare_paper, w=all_atus$weight_adj)  
j<-weighted.mean(all_atus$leisure_paper, w=all_atus$weight_adj) 
k<-weighted.mean(all_atus$tv_paper, w=all_atus$weight_adj) 
l<-weighted.mean(all_atus$socializing_paper, w=all_atus$weight_adj)  
m<-weighted.mean(all_atus$sleeping_paper, w=all_atus$weight_adj)   
n<-weighted.mean(all_atus$ep_paper, w=all_atus$weight_adj)   
o<-weighted.mean(all_atus$otherleisure_paper, w=all_atus$weight_adj)  
p<-weighted.mean(all_atus$other_paper, w=all_atus$weight_adj)   
q<-weighted.mean(all_atus$education_paper, w=all_atus$weight_adj)
r<-weighted.mean(all_atus$civic_paper, w=all_atus$weight_adj) 
s<-weighted.mean(all_atus$ownmedical_paper, w=all_atus$weight_adj) 
d.all.b<- c(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s)

#Then, divide each row by 168-31.325 to get Table 3, column 1.(number of hours a week minus average time spent doing market work)
#Note: do not include the first element (0.231) in table 3.
t3c1<-d.all.b/(168-31.325)

#Columns 2 & 3 (unweighted first differences from 2009 to 2014 in two-year increments)
#Note: below uses slightly different S.E. formula than the original AER paper.
library(plm)
#Bring in the data set:
state_1865_base <- read.dta13("New Data/state_1865_base.dta")
state_1865_base <- state_1865_base[state_1865_base$twoyear>3,] #Here, I'm dropping the data before 2009 (should be 153 obs)

  
out1 <- plm(workapaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear")) #Run the first difference model
se1 <- sqrt(diag(vcovHC(out1, method="arellano", type="HC1", cluster="group"))) #Ask for cluster robust standard errors


out2 <- plm(workupaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear"))
se2 <- sqrt(diag(vcovHC(out2, method="arellano", type="HC1", cluster="group")))
out3 <- plm(childcarepaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear"))
se3 <- sqrt(diag(vcovHC(out3, method="arellano", type="HC1", cluster="group")))
out4 <- plm(homepaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear"))
se4 <- sqrt(diag(vcovHC(out4, method="arellano", type="HC1", cluster="group")))
out5 <- plm(homeproductionpaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear"))
se5 <- sqrt(diag(vcovHC(out5, method="arellano", type="HC1", cluster="group")))
out6 <- plm(homeownpaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear"))
se6 <- sqrt(diag(vcovHC(out6, method="arellano", type="HC1", cluster="group")))
out7 <- plm(shoppingpaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear"))
se7 <- sqrt(diag(vcovHC(out7, method="arellano", type="HC1", cluster="group")))
out8 <- plm(othercarepaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear"))
se8 <- sqrt(diag(vcovHC(out8, method="arellano", type="HC1", cluster="group")))
out9 <- plm(leisurepaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear"))
se9 <- sqrt(diag(vcovHC(out9, method="arellano", type="HC1", cluster="group")))
out10 <- plm(tvpaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear"))
se10 <- sqrt(diag(vcovHC(out10, method="arellano", type="HC1", cluster="group")))
out11 <- plm(socializingpaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear"))
se11 <- sqrt(diag(vcovHC(out11, method="arellano", type="HC1", cluster="group")))
out12 <- plm(sleepingpaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear"))
se12 <- sqrt(diag(vcovHC(out12, method="arellano", type="HC1", cluster="group")))
out13 <- plm(eppaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear"))
se13 <- sqrt(diag(vcovHC(out13, method="arellano", type="HC1", cluster="group")))
out14 <- plm(otherleisurepaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear"))
se14 <- sqrt(diag(vcovHC(out14, method="arellano", type="HC1", cluster="group")))
out15 <- plm(otherpaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear"))
se15 <- sqrt(diag(vcovHC(out15, method="arellano", type="HC1", cluster="group")))
out16 <- plm(educationpaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear"))
se16 <- sqrt(diag(vcovHC(out16, method="arellano", type="HC1", cluster="group")))
out17 <- plm(civicpaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear"))
se17 <- sqrt(diag(vcovHC(out17, method="arellano", type="HC1", cluster="group")))
out18 <- plm(ownmedicalpaper ~ workpaper, data = state_1865_base, model="fd",index=c("idn", "twoyear"))
se18 <- sqrt(diag(vcovHC(out18, method="arellano", type="HC1", cluster="group")))

#Columns 4 & 5 (same as last models, but weighted)
library(multiwayvcov)
dworkpaper<-diff(state_1865_base$workpaper,1)[state_1865_base$twoyear!=6] #pre-difference my main RHS variable
popweight <- state_1865_base$popweight[which(state_1865_base$twoyear!=6)] #extract a vector of proper length/order of weights
idn <- state_1865_base$idn[which(state_1865_base$twoyear!=6)] #extract a vector of proper length/order of cluster ids

dworkapaper<-diff(state_1865_base$workapaper,1)[state_1865_base$twoyear!=6] #pre-difference my LHS variable
out1 <- lm(dworkapaper~ dworkpaper, weights = popweight) #regress the differenced variables, with weights
out1.vcovCL<-cluster.vcov(out1, idn) #cluster the standard errors
coeftest(out1, out1.vcovCL) #Report results. These last four lines of code are repeated 17 times


dworkupaper<-diff(state_1865_base$workupaper,1)[state_1865_base$twoyear!=6]
out2 <- lm(dworkupaper~ dworkpaper, weights = popweight)
out2.vcovCL<-cluster.vcov(out2, idn) 
coeftest(out2, out2.vcovCL)

dchildcarepaper<-diff(state_1865_base$childcarepaper,1)[state_1865_base$twoyear!=6]
out3 <- lm(dchildcarepaper~ dworkpaper, weights = popweight)
out3.vcovCL<-cluster.vcov(out3, idn) 
coeftest(out3, out3.vcovCL)

dhomepaper<-diff(state_1865_base$homepaper,1)[state_1865_base$twoyear!=6]
out4 <- lm(dhomepaper~ dworkpaper, weights = popweight)
out4.vcovCL<-cluster.vcov(out4, idn) 
coeftest(out4, out4.vcovCL)

dhomeproductionpaper<-diff(state_1865_base$homeproductionpaper,1)[state_1865_base$twoyear!=6]
out5 <- lm(dhomeproductionpaper~ dworkpaper, weights = popweight)
out5.vcovCL<-cluster.vcov(out5, idn) 
coeftest(out5, out5.vcovCL)

dhomeownpaper<-diff(state_1865_base$homeownpaper,1)[state_1865_base$twoyear!=6]
out6 <- lm(dhomeownpaper~ dworkpaper, weights = popweight)
out6.vcovCL<-cluster.vcov(out6, idn) 
coeftest(out6, out6.vcovCL)

dshoppingpaper<-diff(state_1865_base$shoppingpaper,1)[state_1865_base$twoyear!=6]
out7 <- lm(shoppingpaper~ dworkpaper, weights = popweight)
out7.vcovCL<-cluster.vcov(out7, idn) 
coeftest(out7, out7.vcovCL)

dothercarepaper<-diff(state_1865_base$othercarepaper,1)[state_1865_base$twoyear!=6]
out8 <- lm(dothercarepaper~ dworkpaper, weights = popweight)
out8.vcovCL<-cluster.vcov(out8, idn) 
coeftest(out8, out8.vcovCL)

dleisurepaper<-diff(state_1865_base$leisurepaper,1)[state_1865_base$twoyear!=6]
out9 <- lm(dleisurepaper~ dworkpaper, weights = popweight)
out9.vcovCL<-cluster.vcov(out9, idn) 
coeftest(out9, out9.vcovCL)

dtvpaper<-diff(state_1865_base$tvpaper,1)[state_1865_base$twoyear!=6]
out10 <- lm(dtvpaper~ dworkpaper, weights = popweight)
out10.vcovCL<-cluster.vcov(out10, idn) 
coeftest(out10, out10.vcovCL)

dsocializingpaper<-diff(state_1865_base$socializingpaper,1)[state_1865_base$twoyear!=6]
out11 <- lm(dsocializingpaper~ dworkpaper, weights = popweight)
out11.vcovCL<-cluster.vcov(out11, idn) 
coeftest(out11, out11.vcovCL)

dsleepingpaper<-diff(state_1865_base$sleepingpaper,1)[state_1865_base$twoyear!=6]
out12 <- lm(dsleepingpaper~ dworkpaper, weights = popweight)
out12.vcovCL<-cluster.vcov(out12, idn) 
coeftest(out12, out12.vcovCL)

deppaper<-diff(state_1865_base$eppaper,1)[state_1865_base$twoyear!=6]
out13 <- lm(deppaper~ dworkpaper, weights = popweight)
out13.vcovCL<-cluster.vcov(out13, idn) 
coeftest(out13, out13.vcovCL)

dotherleisurepaper<-diff(state_1865_base$otherleisurepaper,1)[state_1865_base$twoyear!=6]
out14 <- lm(dotherleisurepaper~ dworkpaper, weights = popweight)
out14.vcovCL<-cluster.vcov(out14, idn) 
coeftest(out14, out14.vcovCL)

dotherpaper<-diff(state_1865_base$otherpaper,1)[state_1865_base$twoyear!=6]
out15 <- lm(dotherpaper~ dworkpaper, weights = popweight)
out15.vcovCL<-cluster.vcov(out15, idn) 
coeftest(out15, out15.vcovCL)

deducationpaper<-diff(state_1865_base$educationpaper,1)[state_1865_base$twoyear!=6]
out16 <- lm(deducationpaper~ dworkpaper, weights = popweight)
out16.vcovCL<-cluster.vcov(out16, idn) 
coeftest(out16, out16.vcovCL)

dcivicpaper<-diff(state_1865_base$civicpaper,1)[state_1865_base$twoyear!=6]
out17 <- lm(dcivicpaper~ dworkpaper, weights = popweight)
out17.vcovCL<-cluster.vcov(out17, idn) 
coeftest(out17, out17.vcovCL)

downmedicalpaper<-diff(state_1865_base$ownmedicalpaper,1)[state_1865_base$twoyear!=6]
out18 <- lm(downmedicalpaper~ dworkpaper, weights = popweight)
out18.vcovCL<-cluster.vcov(out18, idn) 
coeftest(out18, out18.vcovCL)

#Column 6 (same as 4&5, with controls)
dmale<-diff(state_1865_base$male,1)[state_1865_base$twoyear!=6] 
dblack<-diff(state_1865_base$black,1)[state_1865_base$twoyear!=6]
dmarried<-diff(state_1865_base$married,1)[state_1865_base$twoyear!=6] 
dhvchild<-diff(state_1865_base$hvchild,1)[state_1865_base$twoyear!=6] 
dage1<-diff(state_1865_base$age1,1)[state_1865_base$twoyear!=6] 
dage2<-diff(state_1865_base$age2,1)[state_1865_base$twoyear!=6] 
dage3<-diff(state_1865_base$age3,1)[state_1865_base$twoyear!=6] 
dage4<-diff(state_1865_base$age4,1)[state_1865_base$twoyear!=6] 
dage5<-diff(state_1865_base$age5,1)[state_1865_base$twoyear!=6] #omit this due to multicollineariety, as Stata did for published paper
deduc1<-diff(state_1865_base$educ1,1)[state_1865_base$twoyear!=6] 
deduc2<-diff(state_1865_base$educ2,1)[state_1865_base$twoyear!=6] 
deduc3<-diff(state_1865_base$educ3,1)[state_1865_base$twoyear!=6]
deduc4<-diff(state_1865_base$educ4,1)[state_1865_base$twoyear!=6] #omit this due to multicollineariety, as Stata did for published paper


out1 <- lm(dworkapaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight) #regress the differenced variables, with weights
out1.vcovCL<-cluster.vcov(out1, idn) 
coeftest(out1, out1.vcovCL) 

out2 <- lm(dworkupaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight)
out2.vcovCL<-cluster.vcov(out2, idn) 
coeftest(out2, out2.vcovCL)

out3 <- lm(dchildcarepaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight)
out3.vcovCL<-cluster.vcov(out3, idn) 
coeftest(out3, out3.vcovCL)

out4 <- lm(dhomepaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight)
out4.vcovCL<-cluster.vcov(out4, idn) 
coeftest(out4, out4.vcovCL)

out5 <- lm(dhomeproductionpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight)
out5.vcovCL<-cluster.vcov(out5, idn) 
coeftest(out5, out5.vcovCL)

out6 <- lm(dhomeownpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight)
out6.vcovCL<-cluster.vcov(out6, idn) 
coeftest(out6, out6.vcovCL)

out7 <- lm(shoppingpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight)
out7.vcovCL<-cluster.vcov(out7, idn) 
coeftest(out7, out7.vcovCL)

out8 <- lm(dothercarepaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight)
out8.vcovCL<-cluster.vcov(out8, idn) 
coeftest(out8, out8.vcovCL)

out9 <- lm(dleisurepaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight)
out9.vcovCL<-cluster.vcov(out9, idn) 
coeftest(out9, out9.vcovCL)

out10 <- lm(dtvpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight)
out10.vcovCL<-cluster.vcov(out10, idn) 
coeftest(out10, out10.vcovCL)

out11 <- lm(dsocializingpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight)
out11.vcovCL<-cluster.vcov(out11, idn) 
coeftest(out11, out11.vcovCL)

out12 <- lm(dsleepingpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight)
out12.vcovCL<-cluster.vcov(out12, idn) 
coeftest(out12, out12.vcovCL)

out13 <- lm(deppaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight)
out13.vcovCL<-cluster.vcov(out13, idn) 
coeftest(out13, out13.vcovCL)

out14 <- lm(dotherleisurepaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight)
out14.vcovCL<-cluster.vcov(out14, idn) 
coeftest(out14, out14.vcovCL)

out15 <- lm(dotherpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight)
out15.vcovCL<-cluster.vcov(out15, idn) 
coeftest(out15, out15.vcovCL)

out16 <- lm(deducationpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight)
out16.vcovCL<-cluster.vcov(out16, idn) 
coeftest(out16, out16.vcovCL)

out17 <- lm(dcivicpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight)
out17.vcovCL<-cluster.vcov(out17, idn) 
coeftest(out17, out17.vcovCL)

out18 <- lm(downmedicalpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3, weights = popweight)
out18.vcovCL<-cluster.vcov(out18, idn) 
coeftest(out18, out18.vcovCL)

#Column 7
state_1865_base$time1 <- 0                            #Note: Exclude this as reference category
state_1865_base$time1[state_1865_base$twoyear==4] <- 1
state_1865_base$time2 <- 0
state_1865_base$time2[state_1865_base$twoyear==5] <- 1
state_1865_base$time3 <- 0
state_1865_base$time3[state_1865_base$twoyear==6] <- 1

time2 <- state_1865_base$time2[which(state_1865_base$twoyear!=6)] #This code is probably wrong, watch out

#Run the 18 models:
out1 <- lm(dworkapaper~ dworkpaper + time2, weights = popweight) 
out1.vcovCL<-cluster.vcov(out1, idn) 
coeftest(out1, out1.vcovCL) 

out2 <- lm(dworkupaper~ dworkpaper + time2, weights = popweight)
out2.vcovCL<-cluster.vcov(out2, idn) 
coeftest(out2, out2.vcovCL)

out3 <- lm(dchildcarepaper~ dworkpaper + time2, weights = popweight)
out3.vcovCL<-cluster.vcov(out3, idn) 
coeftest(out3, out3.vcovCL)

out4 <- lm(dhomepaper~ dworkpaper + time2, weights = popweight)
out4.vcovCL<-cluster.vcov(out4, idn) 
coeftest(out4, out4.vcovCL)

out5 <- lm(dhomeproductionpaper~ dworkpaper + time2, weights = popweight)
out5.vcovCL<-cluster.vcov(out5, idn) 
coeftest(out5, out5.vcovCL)

out6 <- lm(dhomeownpaper~ dworkpaper + time2, weights = popweight)
out6.vcovCL<-cluster.vcov(out6, idn) 
coeftest(out6, out6.vcovCL)

out7 <- lm(shoppingpaper~ dworkpaper + time2, weights = popweight)
out7.vcovCL<-cluster.vcov(out7, idn) 
coeftest(out7, out7.vcovCL)

out8 <- lm(dothercarepaper~ dworkpaper + time2, weights = popweight)
out8.vcovCL<-cluster.vcov(out8, idn) 
coeftest(out8, out8.vcovCL)

out9 <- lm(dleisurepaper~ dworkpaper + time2, weights = popweight)
out9.vcovCL<-cluster.vcov(out9, idn) 
coeftest(out9, out9.vcovCL)

out10 <- lm(dtvpaper~ dworkpaper + time2, weights = popweight)
out10.vcovCL<-cluster.vcov(out10, idn) 
coeftest(out10, out10.vcovCL)

out11 <- lm(dsocializingpaper~ dworkpaper + time2, weights = popweight)
out11.vcovCL<-cluster.vcov(out11, idn) 
coeftest(out11, out11.vcovCL)

out12 <- lm(dsleepingpaper~ dworkpaper + time2, weights = popweight)
out12.vcovCL<-cluster.vcov(out12, idn) 
coeftest(out12, out12.vcovCL)

out13 <- lm(deppaper~ dworkpaper + time2, weights = popweight)
out13.vcovCL<-cluster.vcov(out13, idn) 
coeftest(out13, out13.vcovCL)

out14 <- lm(dotherleisurepaper~ dworkpaper + time2, weights = popweight)
out14.vcovCL<-cluster.vcov(out14, idn) 
coeftest(out14, out14.vcovCL)

out15 <- lm(dotherpaper~ dworkpaper + time2, weights = popweight)
out15.vcovCL<-cluster.vcov(out15, idn) 
coeftest(out15, out15.vcovCL)

out16 <- lm(deducationpaper~ dworkpaper + time2, weights = popweight)
out16.vcovCL<-cluster.vcov(out16, idn) 
coeftest(out16, out16.vcovCL)

out17 <- lm(dcivicpaper~ dworkpaper + time2, weights = popweight)
out17.vcovCL<-cluster.vcov(out17, idn) 
coeftest(out17, out17.vcovCL)

out18 <- lm(downmedicalpaper~ dworkpaper + time2, weights = popweight)
out18.vcovCL<-cluster.vcov(out18, idn) 
coeftest(out18, out18.vcovCL)

#Column 8:
out1 <- lm(dworkapaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3 + time2, weights = popweight) 
out1.vcovCL<-cluster.vcov(out1, idn) 
coeftest(out1, out1.vcovCL) 

out2 <- lm(dworkupaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3 + time2, weights = popweight)
out2.vcovCL<-cluster.vcov(out2, idn) 
coeftest(out2, out2.vcovCL)

out3 <- lm(dchildcarepaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3 + time2, weights = popweight)
out3.vcovCL<-cluster.vcov(out3, idn) 
coeftest(out3, out3.vcovCL)

out4 <- lm(dhomepaper~ dworkpaper  + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3 + time2, weights = popweight)
out4.vcovCL<-cluster.vcov(out4, idn) 
coeftest(out4, out4.vcovCL)

out5 <- lm(dhomeproductionpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3 + time2, weights = popweight)
out5.vcovCL<-cluster.vcov(out5, idn) 
coeftest(out5, out5.vcovCL)

out6 <- lm(dhomeownpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 + deduc3 + time2, weights = popweight)
out6.vcovCL<-cluster.vcov(out6, idn) 
coeftest(out6, out6.vcovCL)

out7 <- lm(shoppingpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3 + time2, weights = popweight)
out7.vcovCL<-cluster.vcov(out7, idn) 
coeftest(out7, out7.vcovCL)

out8 <- lm(dothercarepaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3 + time2, weights = popweight)
out8.vcovCL<-cluster.vcov(out8, idn) 
coeftest(out8, out8.vcovCL)

out9 <- lm(dleisurepaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3 + time2, weights = popweight)
out9.vcovCL<-cluster.vcov(out9, idn) 
coeftest(out9, out9.vcovCL)

out10 <- lm(dtvpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3 + time2, weights = popweight)
out10.vcovCL<-cluster.vcov(out10, idn) 
coeftest(out10, out10.vcovCL)

out11 <- lm(dsocializingpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3 + time2, weights = popweight)
out11.vcovCL<-cluster.vcov(out11, idn) 
coeftest(out11, out11.vcovCL)

out12 <- lm(dsleepingpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3 + time2, weights = popweight)
out12.vcovCL<-cluster.vcov(out12, idn) 
coeftest(out12, out12.vcovCL)

out13 <- lm(deppaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3 + time2, weights = popweight)
out13.vcovCL<-cluster.vcov(out13, idn) 
coeftest(out13, out13.vcovCL)

out14 <- lm(dotherleisurepaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3 + time2, weights = popweight)
out14.vcovCL<-cluster.vcov(out14, idn) 
coeftest(out14, out14.vcovCL)

out15 <- lm(dotherpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3 + time2, weights = popweight)
out15.vcovCL<-cluster.vcov(out15, idn) 
coeftest(out15, out15.vcovCL)

out16 <- lm(deducationpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3 + time2, weights = popweight)
out16.vcovCL<-cluster.vcov(out16, idn) 
coeftest(out16, out16.vcovCL)

out17 <- lm(dcivicpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3 + time2, weights = popweight)
out17.vcovCL<-cluster.vcov(out17, idn) 
coeftest(out17, out17.vcovCL)

out18 <- lm(downmedicalpaper~ dworkpaper + dmale + dblack + dmarried + dhvchild + dage1 + dage2 + dage3 + dage4 + deduc1 + deduc2 +deduc3 + time2, weights = popweight)
out18.vcovCL<-cluster.vcov(out18, idn) 
coeftest(out18, out18.vcovCL)

